#ifndef _INTERSECT_H_
#define _INTERSECT_H_

#include "../cuCommon/cuTypes.h"
#include "../cuCommon/cuMath.h"

#include "Ray.h"
#include "Box.h"
#include "Triangle.h"
#include "Plane.h"

// -------------------------------------------------------------
// ----------- Ray - Triangle Intersection Code
// -------------------------------------------------------------

SDH_CALL bool Intersect(const STriangle & triangle, 
                        const SRay      & ray, 
                        float3          & tbg)
{
    // intersection point is presumed to be represented through barycentric coordinates:
    // alpha * tri.v0 + beta * tri.v1 + gamma * tri.v2
    // where alpha, beta, gamma >= 0.0f
    // and alpha + beta + gamma = 1
    // so alpha can be calculated as 1 - (beta + gamma)
    // This leads to a linear system with 3 equation and 3 variables

    float3 E1 = triangle.m_V[0] - triangle.m_V[1];  // edge 1
    float3 E2 = triangle.m_V[0] - triangle.m_V[2];  // edge 2
    float3 T  = triangle.m_V[0] - ray.O;

    float p = E2.y * T.z - T.y * E2.z;
    float r = E1.y * T.z - T.y * E1.z;
    float s = E1.y * E2.z - E2.y * E1.z;

    float m = E2.y * ray.D.z - ray.D.y * E2.z;
    float q = ray.D.y * E1.z - E1.y * ray.D.z;

    float inv_denom  = 1.0 / (E1.x * m + E2.x * q + ray.D.x * s);	

    float e3 = E1.x * p - E2.x * r + T.x * s;
    float t = e3 * inv_denom;

    if (t < EPSILON) 
        return (false);

    float n = T.y * ray.D.z - ray.D.y * T.z;

    float e1 = T.x * m - E2.x * n - ray.D.x * p;
    float beta = e1 * inv_denom;

    if (beta < 0.0)
        return (false);

    float e2 = E1.x * n + T.x * q + ray.D.x * r;
    float gamma = e2 * inv_denom;

    if (gamma < 0.0 )
        return (false);

    if (beta + gamma > 1.0)
        return (false);			

    tbg.x = t;
    tbg.y = beta;
    tbg.z = gamma;

    return true;
}

SDH_CALL bool IntersectO(const STriangle & triangle, 
												 const SRay      & ray, 
												 float3          & tbg)
{
	float3 E1 = triangle.m_V[1] - triangle.m_V[0];  // edge 1
	float3 E2 = triangle.m_V[2] - triangle.m_V[0];  // edge 2
	float3 T = ray.O - triangle.m_V[0];
	float3 P = cross(ray.D, E2);
	float3 Q = cross(T, E1);

	float det = dot(P, E1);
	float inv_det = 1.0f / det;

	float t = dot(Q, E2) * inv_det;
	float u = dot(P, T) * inv_det;
	float v = dot(Q, ray.D) * inv_det;

	bool bIntersect = ((t > EPSILON) && (u >= 0.0f ) && (v >= 0.0f) && (u + v <= 1.0f));

	if (bIntersect)
	{
		tbg.x = t;
		tbg.y = u;
		tbg.z = v;
	}

	return bIntersect;
}

SDH_CALL bool IntersectO(const STriangle & triangle, 
												 const SRay      & ray, 
												 float           & tvalue)
{
	float3 E1 = triangle.m_V[1] - triangle.m_V[0];  // edge 1
	float3 E2 = triangle.m_V[2] - triangle.m_V[0];  // edge 2
	float3 T = ray.O - triangle.m_V[0];
	float3 P = cross(ray.D, E2);
	float3 Q = cross(T, E1);

	float det = dot(P, E1);
	float inv_det = 1.0f / det;

	float t = dot(Q, E2) * inv_det;
	float u = dot(P, T) * inv_det;
	float v = dot(Q, ray.D) * inv_det;

	bool bIntersect = ((t > EPSILON) && (u >= 0.0f ) && (v >= 0.0f) && (u + v <= 1.0f));

	if (bIntersect)
	{
		tvalue = t;
	}

	return bIntersect;
}

SDH_CALL bool Intersect(const STriangle & triangle, 
                        const SRay      & ray, 
                        float           & tmin)
{
    // intersection point is presumed to be represented through barycentric coordinates:
    // alpha * tri.v0 + beta * tri.v1 + gamma * tri.v2
    // where alpha, beta, gamma >= 0.0f
    // and alpha + beta + gamma = 1
    // so alpha can be calculated as 1 - (beta + gamma)
    // This leads to a linear system with 3 equation and 3 variables

    float3 E1 = triangle.m_V[0] - triangle.m_V[1];  // edge 1
    float3 E2 = triangle.m_V[0] - triangle.m_V[2];  // edge 2
    float3 T  = triangle.m_V[0] - ray.O;

    float p = E2.y * T.z - T.y * E2.z;
    float r = E1.y * T.z - T.y * E1.z;
    float s = E1.y * E2.z - E2.y * E1.z;

    float m = E2.y * ray.D.z - ray.D.y * E2.z;
    float q = ray.D.y * E1.z - E1.y * ray.D.z;

    float inv_denom  = 1.0 / (E1.x * m + E2.x * q + ray.D.x * s);	

    float e3 = E1.x * p - E2.x * r + T.x * s;
    float t = e3 * inv_denom;

    if (t < EPSILON) 
        return (false);

    float n = T.y * ray.D.z - ray.D.y * T.z;

    float e1 = T.x * m - E2.x * n - ray.D.x * p;
    float beta = e1 * inv_denom;

    if (beta < 0.0)
        return (false);

    float e2 = E1.x * n + T.x * q + ray.D.x * r;
    float gamma = e2 * inv_denom;

    if (gamma < 0.0 )
        return (false);

    if (beta + gamma > 1.0)
        return (false);			

    tmin = t;

    return true;
}

SDH_CALL bool Intersect(const float3 & v0, 
												const float3 & v1, 
												const float3 & v2,
                        const float3 & orig, 
												const float3 & dir,
                        float        & tmin)
{
	//@ intersection point is presumed to be represented through barycentric coordinates:
	//@ alpha * tri.v0 + beta * tri.v1 + gamma * tri.v2
	//@ where alpha, beta, gamma >= 0.0f
	//@ and alpha + beta + gamma = 1
	//@ so alpha can be calculated as 1 - (beta + gamma)
	//@ This leads to a linear system with 3 equation and 3 variables

	float3 E1 = v1 - v0;  // edge 1
	float3 E2 = v2 - v0;  // edge 2
	float3 P = cross(dir, E2);
	float det = dot(P, E1);

	float3 T = orig - v0;
	float3 Q = cross(T, E1);

	float inv_det = 1.0f / det;

	float t = dot(Q,  E2) * inv_det;
	if (t <= EPSILON) return false;

	float u = dot(P,   T) * inv_det;
	float v = dot(Q, dir) * inv_det;

	bool bIntersect = ((u >= 0.0f ) && (v >= 0.0f) && (u + v <= 1.0f));

	if (bIntersect)
	{
		tmin = t;
	}

	return bIntersect;
}

SDH_CALL bool Intersect(const STriangle & triangle, 
                        const float3    & point, 
                        float3          & tbg)
{
    SPlane plane;

    float3 E1 = triangle.m_V[1] - triangle.m_V[0];  // edge 1
    float3 E2 = triangle.m_V[2] - triangle.m_V[0];  // edge 2

    float3 triangleNormal = cross(E1, E2);
    Plane(&plane, triangleNormal, triangle.m_V[0]);

    float3 planePoint;
    float t = SignedDistanceToPoint(plane, point, planePoint);

    float3 v2 = planePoint - triangle.m_V[0]; 

    float E1E1 = dot(E1,E1);
    float E1E2 = dot(E1,E2);
    float E2E2 = dot(E2,E2);

    float E1v2 = dot(E1,v2);
    float E2v2 = dot(E2,v2);

    float invDenom = 1.0f / (E1E1 * E2E2 - E1E2 * E1E2);
    float b = (E2E2 * E1v2 - E1E2 * E2v2) * invDenom;
    float g = (E1E1 * E2v2 - E1E2 * E1v2) * invDenom;
   
    bool bIntersect = ((b >= 0.0f ) && (g >= 0.0f) && (b + g <= 1.0f));

    if (bIntersect)
    {
        tbg.x =  abs(t);
        tbg.y =  b;
        tbg.z =  g;
    }
			
    return bIntersect;
}

SDH_CALL bool Intersect(const STriangle & triangle, 
                        const float3    & point, 
                        float           & tmin)
{
    SPlane plane;

    float3 E1 = triangle.m_V[1] - triangle.m_V[0];  // edge 1
    float3 E2 = triangle.m_V[2] - triangle.m_V[0];  // edge 2

    float3 triangleNormal = cross(E1, E2);
    Plane(&plane, triangleNormal, triangle.m_V[0]);

    float3 planePoint;
    float t = SignedDistanceToPoint(plane, point, planePoint);

    float3 v2 = planePoint - triangle.m_V[0]; 

    float E1E1 = dot(E1,E1);
    float E1E2 = dot(E1,E2);
    float E2E2 = dot(E2,E2);

    float E1v2 = dot(E1,v2);
    float E2v2 = dot(E2,v2);

    float invDenom = 1.0f / (E1E1 * E2E2 - E1E2 * E1E2);
    float b = (E2E2 * E1v2 - E1E2 * E2v2) * invDenom;
    float g = (E1E1 * E2v2 - E1E2 * E1v2) * invDenom;

    bool bIntersect = ((b >= 0.0f ) && (g >= 0.0f) && (b + g <= 1.0f));

    if (bIntersect)
    {
        tmin =  abs(t);
    }

    return bIntersect;
}

// -------------------------------------------------------------
// ----------- Ray - Box Intersection Code
// -------------------------------------------------------------

SDH_CALL 
bool Intersect(const S_Box & box, const SRay & ray)
{
    float txmin, txmax, tymin, tymax, tzmin, tzmax;

    float inv_dir_x = 1.0f / ray.D.x;
    if (ray.D.x >= 0) 
    {
        txmin = (box.pmin.x - ray.O.x) * inv_dir_x;
        txmax = (box.pmax.x - ray.O.x) * inv_dir_x;
    }
    else 
    {
        txmin = (box.pmax.x - ray.O.x) * inv_dir_x;
        txmax = (box.pmin.x - ray.O.x) * inv_dir_x;
    }

    float inv_dir_y = 1.0f / ray.D.y;
    if (ray.D.y >= 0) 
    {
        tymin = (box.pmin.y - ray.O.y) * inv_dir_y;
        tymax = (box.pmax.y - ray.O.y) * inv_dir_y;
    }
    else 
    {
        tymin = (box.pmax.y - ray.O.y) * inv_dir_y;
        tymax = (box.pmin.y - ray.O.y) * inv_dir_y;
    }

    float inv_dir_z = 1.0f / ray.D.z;
    if (ray.D.z >= 0) 
    {
        tzmin = (box.pmin.z - ray.O.z) * inv_dir_z;
        tzmax = (box.pmax.z - ray.O.z) * inv_dir_z;
    }
    else 
    {
        tzmin = (box.pmax.z - ray.O.z) * inv_dir_z;
        tzmax = (box.pmin.z - ray.O.z) * inv_dir_z;
    }

    float t0, t1;

    t0 = (txmin > tymin) ? txmin: tymin;
    t0 = (tzmin > t0) ? tzmin: t0;

    t1 = (txmax < tymax) ? txmax: tymax;
    t1 = (tzmax < t0) ? tzmax: t1;

    return ( (t0 < t1) && (t1 > EPSILON) );
}

SDH_CALL 
bool Intersect(const S_Box & box, const SRay & ray, float & t0, float & t1)
{
    float txmin, txmax, tymin, tymax, tzmin, tzmax;

    float inv_dir_x = 1.0f / ray.D.x;
    if (ray.D.x >= 0) 
    {
        txmin = (box.pmin.x - ray.O.x) * inv_dir_x;
        txmax = (box.pmax.x - ray.O.x) * inv_dir_x;
    }
    else 
    {
        txmin = (box.pmax.x - ray.O.x) * inv_dir_x;
        txmax = (box.pmin.x - ray.O.x) * inv_dir_x;
    }

    float inv_dir_y = 1.0f / ray.D.y;
    if (ray.D.y >= 0) 
    {
        tymin = (box.pmin.y - ray.O.y) * inv_dir_y;
        tymax = (box.pmax.y - ray.O.y) * inv_dir_y;
    }
    else 
    {
        tymin = (box.pmax.y - ray.O.y) * inv_dir_y;
        tymax = (box.pmin.y - ray.O.y) * inv_dir_y;
    }

    float inv_dir_z = 1.0f / ray.D.z;
    if (ray.D.z >= 0) 
    {
        tzmin = (box.pmin.z - ray.O.z) * inv_dir_z;
        tzmax = (box.pmax.z - ray.O.z) * inv_dir_z;
    }
    else 
    {
        tzmin = (box.pmax.z - ray.O.z) * inv_dir_z;
        tzmax = (box.pmin.z - ray.O.z) * inv_dir_z;
    }

    t0 = (txmin > tymin) ? txmin: tymin;
    t0 = (tzmin > t0) ? tzmin: t0;

    t1 = (txmax < tymax) ? txmax: tymax;
    t1 = (tzmax < t0) ? tzmax: t1;

    return ( (t0 < t1) && (t1 > EPSILON) );
}

SDH_CALL 
bool Intersect(const float3 & pmin,
							 const float3 & pmax,
							 const float3 & rayO,
							 const float3 & rayD,
							 float & t0, 
							 float & t1)
{
    float txmin, txmax, tymin, tymax, tzmin, tzmax;

    float inv_dir_x = 1.0f / rayD.x;
    if (rayD.x >= 0) 
    {
        txmin = (pmin.x - rayO.x) * inv_dir_x;
        txmax = (pmax.x - rayO.x) * inv_dir_x;
    }
    else 
    {
        txmin = (pmax.x - rayO.x) * inv_dir_x;
        txmax = (pmin.x - rayO.x) * inv_dir_x;
    }

    float inv_dir_y = 1.0f / rayD.y;
    if (rayD.y >= 0) 
    {
        tymin = (pmin.y - rayO.y) * inv_dir_y;
        tymax = (pmax.y - rayO.y) * inv_dir_y;
    }
    else 
    {
        tymin = (pmax.y - rayO.y) * inv_dir_y;
        tymax = (pmin.y - rayO.y) * inv_dir_y;
    }

    float inv_dir_z = 1.0f / rayD.z;
    if (rayD.z >= 0) 
    {
        tzmin = (pmin.z - rayO.z) * inv_dir_z;
        tzmax = (pmax.z - rayO.z) * inv_dir_z;
    }
    else 
    {
        tzmin = (pmax.z - rayO.z) * inv_dir_z;
        tzmax = (pmin.z - rayO.z) * inv_dir_z;
    }

    t0 = (txmin > tymin) ? txmin: tymin;
    t0 = (tzmin > t0) ? tzmin: t0;

    t1 = (txmax < tymax) ? txmax: tymax;
    t1 = (tzmax < t0) ? tzmax: t1;

    return ( (t0 < t1) && (t1 > EPSILON) );
}
#endif
